%% Script for analyzing the Structured Illumination based FTH experiments 
% Double-Constraint Gerchberg-Saxton CGH Synthesis;
% 
% Written by Kahraman Keskinbora kahraman@mit.edu keskinbora@is.mpg.de
% in 2021 for the paper "Maskless Fourier Transform Holography"
%
%
%
%% Initialize the workspace

clear all
close all
clc

load cmap.mat % this is a color scheme used in some of the plots

%% Units;

meters      =1;
microns     =1e-6 * meters;
nanometers  =1e-9 * meters;
cm          =10^-2 * meters;

seconds     =1;
miliseconds =1e-6 * seconds;
nanoseconds =1e-9 * seconds;

hertz       =1/seconds;
kilohertz   =1e6 * hertz;
gigahertz   =1e9 * hertz;
terahertz   =1e12* hertz;
petahertz   =1e15* hertz;

eV          =1;
keV         =eV * 1000;
meV         =eV / 1000;


%% Constants;

c0          = 299792458 * meters/seconds;
e0          = 8.854187817e-12 * 1/meters;
u0          = 1.2566370614e-6 * 1/meters;

%% Load data
% Please correct this path for your computer
fid = imread('D:\PostDoc\MIT_DFG_Work\StIXH_Experiments\Table_Top\Experimental_Data\071921_qpo_SS1_1.tiff');
data = fid;

[Xsize,Ysize]=size(fid);


%% Create a reference for affine transform, crop and center the whole deck to remove the phase ramp

bg_ref = imref2d(size(data));


shiftX = 25;            % Distance of the brightest spot from 512 in x in pix
shiftY = -30;           % Distance of the brigthest spot from 512 in y in pix

T = [1 0 0; 0 1 0; shiftX shiftY 1];
tform = affine2d(T);


slice = data;
[shifted_slice,shifted_ref] = imwarp(data,tform,'OutPutView',bg_ref);
shifted_data = shifted_slice(:,[129:1152]);

   
%% Reconstruct

recons_data = fftshift(fft2(fftshift(shifted_data)));
recons = recons_data;

% Calculate the coordinates 

lambda  = 635*nanometers;                   % Wavelenght

Det_z = 5*cm;                               % Sample to detector distance
N     = Xsize;                              % Data size
dp_det= 5.2*microns;                        % Detector pixel size
w     = c0/lambda;                          % Frequency
dr    = (2*pi*Det_z*c0)/(N*dp_det*w);       % Sampling (pixel) size in the sample plane
L2    = 1024*dp_det;                        % Detector plane side length
d2    = -L2/2:dp_det:L2/2-dp_det;           % Detector coordinates

dx1   = (lambda*Det_z)/L2;                  % Sample coordinate pixel size
L1    = 1024*dx1;                           % Sample side length

d1    = -L1/2:dx1:L1/2-dx1;                 % Sample coordinates

%% Plot results

amp_image = abs(recons([559:633],[559:633]));
phase_image = angle(recons([559:633],[559:633]));

N_crop = size(amp_image);
dcrop = 0:dx1:N_crop*dx1-dx1;

figure(1)
imagesc(dcrop*1e6,dcrop*1e6,sqrt(amp_image))
axis square
axis tight
% xlim([250 750])
% ylim([250 750])
xlabel 'Distance \mum'
ylabel 'Distance \mum'
colormap 'bone'
% caxis([6.3 7.3])
colorbar 

figure(2)
imagesc(dcrop*1e6,dcrop*1e6,phase_image)
axis square
axis tight
% xlim([250 750])
% ylim([250 750])
xlabel 'Distance \mum'
ylabel 'Distance \mum'
colormap(cmap)
colorbar 

figure(3)
imagesc(d2*1e3,d2*1e3,log10(shifted_data))
axis square
axis tight
xlim([-2.5 2.5])
ylim([-2.5 2.5])
xlabel 'Distance mm'
ylabel 'Distance mm'
colormap 'jet'
colorbar 

% Calculate the line profile

d_profile = 0:dx1:20*dx1-dx1;
lineprofile1 = phase_image((1:20),56);
lineprofile2 = phase_image((1:20),57);
line_average = (lineprofile1+lineprofile2)/2;
figure(4)
plot(d_profile*1e6,lineprofile1,'bo',d_profile*1e6,lineprofile2,'g*',d_profile*1e6,line_average,'k','LineWidth',2)
legend
ylabel '\Phi radians'
xlabel 'Distance pixels'
legend('Line Profile 1', 'Line Profile 2', 'Average')

delta_phi = (max(line_average)-min(line_average))

